Thermocapillary migration of a planar droplet 
at moderate and large Marangoni numbers 



Zuo-Bing Wu*^ and Wen-Rui Hu+ 
State Key Laboratory of Nonlinear Mechanics* 
and National Microgravity Laboratory +, 
Institute of Mechanics, Chinese Academy of Sciences, 
Beijing 100190, China 

December 14, 2011 



1 Corresponding author. Tel:. +86-10-82543955; fax.: -^86-10-82543977. Email address: 
wuzb@lnni.iinech.ac.cn (Z.-B. Wu). 

1 



Abstract 

Thermocapillary migration of a planar non-deformable droplet in 
flow fields with two uniform temperature gradients at moderate and 
large Marangoni numbers is studied numerically by using the front- 
tracking method. It is observed that the thermocapillary motion of 
planar droplets in the uniform temperature gradients is steady at mod- 
erate Marangoni numbers, but unsteady at large Marangoni numbers. 
The instantaneous migration velocity at a fixed migration distance 
decreases with increasing Marangoni numbers. The simulation results 
of the thermocapillary droplet migration at large Marangoni num- 
bers are found in qualitative agreement with those of experimental 
investigations. Moreover, the results concerned with steady and un- 
steady migration processes are further confirmed by comparing the 
variations of temperature fields inside and outside the droplet. It 
is evident that at large Marangoni numbers the weak transport of 
thermal energy from outside of the droplet into inside cannot satisfy 
the condition of steady migration process, which implies that the ad- 
vection around the droplet is a more significant mechanism for heat 
transfer across/around the droplet at large Ma numbers. Furthermore, 
from the condition of overall steady-state energy balance in the flow 
domain, the thermal flux across its surface is studied for a steady ther- 
mocapillary droplet migration in a flow field with uniform temperature 
gradient. By using the asymptotic expansion method, a nonconserva- 
tive integral thermal flux across the surface is identified in the steady 
thermocapillary droplet migration at large Marangoni numbers. This 
nonconservative flux may well result from the invalid assumption of 
quasi-steady state, which indicates that the thermocapillary droplet 
migration at large Marangoni numbers cannot reach steady state and 
is thus a unsteady process. 
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1 Introduction 



The transport phenomenon of droplets/bubbles in a hquid is a very im- 
portant topic for both fundamental hydrodynamics and practical applica- 
tions such as production of pure materials in manufacturing industry and 
mass transfer in chemical engineering. Under normal gravity, the motion of 
droplet /bubbles results from the buoyancy when the densities of two fluids 
are different. With fast development of space exploration, the studies on 
the physical mechanism of droplet/bubble migration phenomena under re- 
duced gravity become more and more important. In this case, the buoyant 
effect vanishes, the droplet/bubble moves as a result of the variance of in- 
terface tension. Thus, in the microgravity environment, a droplet/bubble 
suspended in an ambient fluid will move in the direction of temperature 
gradient due to thermocapillary force [1]. Thermocapillary motion of a single 
droplet was firstly examined both theoretically and experimentally by Young, 
Goldstein & Block (1959) [2J. They gave an analytical prediction on its mi- 
gration speed in the limit case of zero Reynolds (Re) and zero Marangoni 
(Ma) numbers, which is called as YGB model. Since then, the thermo- 
capillary migration of a bubble has been studied extensively by a series of 
theoretical analyses[3l HI El |6], numerical simulations [71 El |9l [10] and experi- 
mental investigations [1 1 j . In the mean time, several numerical techniques for 
treating the two-phase flow, such as the front-tracking method [T^ and 
the level-set method[ll], have also been developed, which may provide effec- 
tive techniques to directly investigate thermocapillary migration processes 
of bubbles or droplets [151 [IHl [HI [18], interfacial mass transfer [T9l [20] and 
interfacial flows with soluble surf act ant s [^ [2^ . 

For the migration of a droplet, the experimental result for the migration 
speed at small Re numbers obtained by Braun et al(1993)[23] agrees with the 
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YGB model. To include small inertial effects, the YGB model analysis was 
extended to the range of small Ma numbers[24j. For finite Ma numbers, sev- 
eral numerical simulation on the three-dimensional thermocapillary motion of 
non-deformable and deformable droplets were reported by Wang et al[25] and 
Haj-Hariri et al[T7], respectively. They used the front-tracking and level-set 
methods to catch the interface and investigate the effects of physical param- 
eters on migration speeds and mobility, respectively. For large Ma numbers, 
Balasubramanian & Subramanian(2000)[2B] used thermal boundary layers 
and found that the migration speed of a droplet increased with increasing 
Ma number, which is in qualitative agreement with corresponding numerical 
simulations |27]. Both the theoretical analysis and numerical simulation are 
based on the assumptions of quasi-steady state and non-deformation of the 
droplet. However, the experimental results of Hadland et al(1999)[TT] and 
Xie et al(2005)[28] were not in qualitative agreement with the above theo- 
retical and numerical results, and it was shown that the droplet migration 
speed non-dimensionalized by the YGB velocity decreased as Ma number 
increased. The experiment investigation was completed in several ranges of 
large Ma numbers, where the droplet migration was in an accelerating state 
and did not reach a steady one. Recently, a numerical investigation based on 
an axisymmetric droplet model[29l|30] found that the steady thermocapillary 
migration process did exist in a laboratory coordinate system, and verified 
the above experimental results for the case of large Ma numbers, however, 
the effects of Capillary numbers were not given in the calculations. Owing 
to the invariance theory under transformation between two inertia frames, 
the above numerical result in a laboratory coordinate frame, i.e., the steady 
migration speed of the droplet decreases as Ma increases, should be in agree- 
ment with the theoretical and numerical results in a reference frame moving 
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with the steady droplet velocity EZ]- However, this seems impossible. 
Moreover, Herrmann et al[TB] adopted numerical method to investigate the 
thermocapillary motion of deformable droplets and indicated that for large 
Ma numbers the assumption of quasi-steady state was not valid. Therefore, 
the thermocapillary droplet migration at large Ma numbers is still a topic to 
be studied with emphasis laid on its physical mechanism. 

The planar or cylindrical droplet/bubble as a simple model has been 
extensively used to study its dynamical mechanism [211 [321 E3]- In this paper, 
we use the front-tracking method to numerically study the thermocapillary 
migration of a non- deformable planar droplet in the liquids at moderate and 
large Ma numbers, and analyze in detail the relation of the migration velocity 
to the temperature distribution inside and outside the droplet. Moreover, by 
using asymptotic analysis, we investigate the continuity of thermal flux across 
the surface based on the overall energy balance of the droplet, and analyze 
the existence of quasi-steady migration of the droplet at large Ma numbers. 

2 Governing equations 

Consider the thermocapillary migration of a planar droplet in a continu- 
ous phase fluid of infinite extent under a uniform temperature gradient G. 
Gravity and deformation of the droplet shape are ignored. Two-dimensional 
continuous, momentum and energy equations for the continuous phase fluid 
and the droplet in a laboratory coordinate system are written as follows 

^ + V ■ (p.v,) = 0, 

^ + V ■ (p.v.Vi) = -Vpi + V ■ /i.(Vv, + V^Vi + F,, (1) 
^ + V-(v,T,) = |V-(A;iVTO, 
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where Vj and Tj are velocity and temperature, respectively. F^- is the surface 
tension force acting on the interface, pi, fii, ki, Ki represent density, dynamic 
viscosity, thermal conductivity, and thermal diffusivity, respectively Sym- 
bols with subscript 1 and 2 denote physical coefficients of the continuous 
fluid and the droplet, respectively. The solutions of Eqs. (1) have to satisfy 
the boundary conditions at infinity 

vi = 0,ri ^To + G^, (2) 

where Tq is the undisturbed temperature of the continuous phase and the 
boundary conditions at the interface (r^, z^^ of the two fiuids 

Ti{n,Zh,t)^T2{n,Zh,t), (3) 

where n is a unit vector normal to the interface. In the modelling assump- 
tions, both fiuids are immiscible and the physical properties are constant. 
The equations of state for density, viscosity, heat conduction and heat diffu- 
sivity are written as follows 

d£i^ _ djM _ dh _ _ Q 
dt dt dt dt 

The reference velocity is defined as 

-arGRo/fJ.!, (5) 

where Rq is the radius of the droplet, and ari— da/dT) is the change rate 
of interfacial tension with temperature. By taking Rq, Vo and GRq as the 
characteristic quantities to make coordinates, velocity and temperature di- 
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mensionless, Eqs. (1) are rewritten in the following non-dimensional form 
V • = 0, 

+ V • (av,v,) = -Vp, + • /x,(Vv, + V^v,) + f,, (6) 
f + V-(v,T,) = ^|V-(/c,VT,), 

where the physical coefficients are non-dimensionlized by the characteristic 
quantities of continuous fluid and — F^Rq/ piVg. The boundary conditions 
(2) and (3) are non-dimensionlized as 

vi = o,ri^ro + z, (7) 

at infinity and 

yi(rb,Zb,t) = V2(rb,Zb,t), 

Ti {n, Zb,t) = T2 {rb, Zb,t), (8) 

h^{rb,Zb,t) = k2^{rb,Zb,t), 
at the interface between two fluids. The Reynolds number and Marangoni 
number are defined as 

Re^P^, Ma^'-^^PrRe, (9) 

where Pr = Hi/piKi is the Prandtl (Pr) number. In what follows, the undis- 
turbed temperature Tq and non-dimensional physical parameters (pi — /ii — 
ki — Ki — 1) of the continuous phase are reduced for simplicity, except when 
otherwise indicated. 
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3 Numerical simulation of thermocapillary droplet 
migration at moderate and large Ma num- 
bers 

3.1 Models and numerical methods 

As shown schematically in Fig.l, the symmetric axis of the container is taken 
as the z-axis. A droplet is placed initially at the center of coordinates and 
then moved along the 2;-axis. Thus, the solution of Eqs. (6) should satisfy the 
following initial conditions in the whole domain x e [xq, xi] and z e [zq, zi] 

V, = 0, Ti^z (10) 

and non-slip/periodic boundary conditions at the top and bottom walls/the 
horizontal boundaries 

Vi{x,Zo) =y^iix,zi) =0, Ti{x,Zo) = zo,Ti{x,Zi) = zi, ^^^^ 
Vi{xo, z) = vi(a;i, z), Ti{xq, z) = Ti{xi, z). 

In the computation, we use a fixed regular staggered MAC grid in the 
computational domain. To discretize Eqs. (6), we adopt a second-order 
central difference scheme for the spatial variables and an explicit predictor- 
corrector second-order scheme for time integration. The predictor-corrector 
method is a combination of the explicit Euler and the implicit trapezoidal 
methods to obtain an explicit technique with better convergence character- 
istic. In the method, the solution at time step n -|- 1 is predicted by using 
the explicit Euler method 

0Ul=r + /(^n,r)At, (12) 

where f indicates that this is not the final value of the solution at t„+i. 
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Rather, the solution is corrected by applying the trapezoid rule 

r^' = r + lifituA) + f{tn+iAi+i)W. (13) 

To achieve the second-order accurate time integration of the velocity and 
temperature fields in Eqs. (6), we employ the Chorin's projection method to 
outline the first-order Euler integration in (12) as follows. 

Since both fiuids are assumed immiscible, all physical coefficients are 
discontinuous across the interface. The interface is captured and updated 
by the front-tracking method IT3]. When the interface is moved to a new 
position, the density field is updated. The interface is considered to have a 
finite width so that the density across the interface is continuous. Here, a 
weighting function suggested by Peskin[M] is adopted as 



>p) = d{xp - iAx)d{zp - jAz), (14) 



where 



(l/4Ar)[l + cos(7rr/2/i)], |r| < 2Ar, 
d{r) = { (15) 
0, |r| > 2Ar, 

and (xp, Zp) is the interface node. Once the density is updated, the velocity 

field will be computed by the Chorin's projection method, which is divided 

into two steps. One is a prediction step, where the effect of the pressure is 

ignored 

- — = -V ■ (prv>r) + ■ /ir(vvr + v^v^ + f.. (le) 

Another is a correction step in terms of the pressure gradient 

n+1 n+1 _ ^n+1 * 

- — '-^^ — - = -vpr\ (17) 
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where the pressure is obtained by solving the following Poisson equation 

V ■ -irVpr' = • V*. (18) 

In solving Eq. (18), we use the successive over relaxation iteration method 
to get p^^^. When the pressure is obtained, the corrected velocity field v""*"^ 
is determined from Eq. (17). Similarly, the energy equation is discretized in 
the form 

rpn+l rpn -| n 

^— i- = -V ■ (v,"+i7;") + — T^V ■ iM^Tt). (19) 

In terms of the corrected velocity field v"'*'^, the temperature field 77*+^ is 
determined. Until now, by using the projection method, the first-order time 
integration (j)^n+i of the velocity, pressure and temperature fields is completed. 
In the mean time, other physical coefficients (//, /c, k) across the interface at 
the time step n + 1 are also updated to have the same distribution as the 
density. Repeating the above process, we get a second first-order accurate 
solution <i)'n+2 time step n-\- 2 based on the first-order accurate solution 
+1 

0].V2 = + fUn+u 0I+1) Ai. (20) 



Finally, the solution for the second-order time integration is obtained as 
follows 

= (0- + 0tt^,)/2. (21) 

Since the droplet in the migration process is assumed non-deformable, the 
vertical area average velocity in the droplet is taken as the droplet migration 
velocity V^. The nodes of interface are moved with this velocity at each time 
step. 

In solving Eq. (16), the surface tension force fo- is determined by referring 
to the temperature field. In general, the surface tension force on a short front 
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element is defined as 



J As 



I — (at)ds=((7t)2-((7t)i=A2i((7t), 



(22) 



where t is an unit tangent vector, s is the arc length along the interface, a 
is the surface tension coefficient written as 



where is the surface tension coefficient at a reference temperature To, and 
Gt is a negative constant for most fluids. By using the above characteristic 
quantities, the non-dimensional surface tension force is written in the form 
of body force as 



where C&{=vqi1i/ uq) is the Capillary number. To calculate the surface ten- 
sion force fo-, the surface temperature on the interface is firstly obtained by 
interpolating values on the grid points. The tangent vector is computed from 
a Lagrangian polynomial fitting through four interface nodes. Then, the sur- 
face tension force on the interface is distributed to the grid points by means 
of weighting function (14). 

3.2 Results and analysis 

To check the sensitivity of the results to grid refinements, we perform calcula- 
tions for a planar droplet migration at Re=5, Ma=20, Ca=0. 01666 and p2— 
H2— k2— «;2=0.5 using the method described above. The computational do- 
main is chosen as 4 x 8. Based on 64 x 128, 96 x 192 and 128 x 256 grid points. 



(23) 



6f, = 5¥,{R^/pivl)/{Rl5x5z) 
= ^2i{cr^) / PivIRq5x5z 
= A2i[((7o/'i;o/ii - T)t\/Re5x5z 
= A2i[(l/Ca - T)t\lRe5x5z, 



(24) 
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i.e., 16, 24 and 32 grid points per droplet radius, the time evolution of the 
droplet migration velocity is calculated and plotted in Fig. 2. The migration 
velocity curve seems to converge when the grid becomes finer. The difference 
in the migration velocities computed with 24 and 32 grid points per droplet 
radius is very small (about 1.5%). In the following calculations we fix 24 grid 
points per droplet radius as the grid resolution. To further validate our code, 
we compare the current computation results with Nas & Tryggvason's|15]. 
where the deformation of the droplet is considered. In Fig. 3, it is observed 
that both results have the same trends and the migration velocities are close 
together. 

In the following calculations, we adopt the silicone oil of nominal viscosity 
5cst and the FC-75 Fluorinert liquid, i.e., the working media in the space 
experiment [28] , as the continuous phase fiuid and the droplet, respectively. 
The physical parameters of the continuous fiuid and the droplet at temper- 
ature 25°C are given in Table I. ax is fixed as -0.044 dyn / cmK and o"o ~ 
6 dyn/cm|35] is adopted. From the values of the continuous fiuid parame- 
ters, the Pr number and the capillary length Aq = \J (^q/ PiQo (with the earth's 
gravity = 980cm/ s"^) are determined as 67.8 and 0.08cm, respectively. The 
most of droplets in the space experiment |28j have Rq > Aq, which refers to 
the domination of the gravitational effect to the droplet shapes on the earth. 
However, in the microgravity environment (the effective gravity Qe is about 
0(10^^) of go), the gravitational effect is neglected(-Ro -^e), so the droplet 
shapes are dominated by the capillary effect. The computational domain is 
chosen as {x, z} G {[—4, 4], [—4, 16]} and the resolution is fixed at 192 x 480. 
The initial droplet is placed at the position (0,0) and the time step is 0.0002. 
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3.2.1 Flow field with the temperature gradient G=12 K/cm 

In the space experiments with G=12 K/cmp8], Re and Ma lie respectively 
in the ranges of 4.5-302.6 and 145-5525, their specific values depending on 
Rq. To simulate the experimental processes, the physical coefficients in the 
droplet migration processes are determined by changing Rq. The correspon- 
dence of Re, Ma and Ca to Rq is presented in Table II, where Re is in a range 
of moderate values and Ma have both moderate and large numbers. Fig. 4 
displays the time evolution of droplet migration velocities for five sets of non- 
dimensional coefficients. In the present range of Ma, the migration velocities 
versus time have complex behaviors, which can be classified into three types 
based on the curve characters. At Ma=44.7(Re=0.66), the initial migration 
velocity increases sharply before t = 3, and then drops to approach a steady 
value. For Ma=402. 5-1118. l(Re=5. 93-16. 5), the initial accelerating process 
has smaller peak value as Ma increases. After the increasing-decreasing os- 
cillation process, the terminal droplet migration velocity increases with time, 
i.e., the droplet migration is in an accelerating state. The slope of the curve 
increases as Ma increases. For Ma=2191. 6-3622. 8(Re=32. 3-53. 4), the droplet 
migration velocity increases monotonously with time and decreases with in- 
creasing Ma. We can thus conclude that in the time frame under investigation 
the thermocapillary droplet migration is steady at moderate Ma numbers, 
but becomes unsteady at large Ma numbers. In the two space experiments. 
Figs. 4 of [11] and [2E] showed that the whole migration processes were un- 
steady and didn't reach any steady state. Even a plateau appears in the curve 
of migration velocity vs migration distance, the migration process seems to 
be an accelerating one after the slow varying period. To further compare 
with the experimental investigation [28] . we take several fixed migration dis- 
tances Iz and determine the relation between instantaneous non-dimensional 
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migration velocity Vz and Ma numbers. The numerical and experimental re- 
sults are plotted in Fig. 5, from which it is evident that both the numerical 
and experimental migration velocities of droplet decrease as Ma increases 
in the range of large Ma numbers. Hence, at large Ma numbers, the above 
simulation results are in qualitative agreement with those of experimental 
investigations [m 128] . 

In order to understand the phenomena exhibited in droplet migration 
processes, it is important to analyze the evolution of the velocity and tem- 
perature fields. Fig. 6 displays the computed velocity fields at t = 20 in 
both the laboratory coordinate frame and the reference frame moving with 
the droplet at Re=16.5(Ma=1118.1). In the laboratory coordinate frame, 
the streamlines for a moving droplet are closed and symmetric about the 
2;-axis. In the reference frame, when the external streamlines go around the 
droplet, a pair of vortices is formed inside the droplet. It is evident that in 
the reference frame recirculation flows in both the continuous phase fluid and 
the droplet are driven by the surface tension force generated by the temper- 
ature gradient along the surface. In Fig. 7, we depict the pattern evolution 
of streamlines with time in a reference frame moving with the droplet at 
Re=16.5(Ma=1118.1). Initially, there appear two vortices symmetric about 
the vertical diameter. Along with the rising of the droplet, the pair of vor- 
tices in the droplet is kept, but the vortex centers are moving up. In the 
whole process, the basic types of streamlines for both the internal and exter- 
nal motions are not changed. For moderate Re numbers, both the convection 
and viscous terms in the momentum equation have important effects on the 
fluid flow. The external flow just passes around the droplet and does not 
separate from the droplet surface, and thus the computed velocity fields are 
similar in the range of moderate Re numbers. 
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In the range of Rq, the fluid flow has behaviors similar to those for mod- 
erate Re numbers, but the thermal transfer exhibits different characters for 
moderate and large Ma numbers. For moderate Ma numbers, both the heat 
convection and the heat conduction have important effects on the energy 
transfer. Fig. 8 displays the time evolution of isotherms at moderate Ma 
(Ma=44.7, Re=0.66), which corresponds to that of migration velocities shown 
in Fig. 4. As given in Eq. (23), the surface tension coefficient decreases with 
the increasing of the local temperature. For a temperature field with its 
gradient in the z direction, the generated surface tension force is a net force 
along the surface. At the beginning, the droplet starts to move towards the 
warm side under the action of net force. It induces in turn viscous stresses 
in both fluids, which causes streamlines inside and outside the droplet to 
form double vortices and to go around the droplet, respectively. The tem- 
perature held inside the droplet is affected by the two rotating vortices. The 
horizontal isotherm T = in the droplet is moving up, as well as bending 
along the migration direction. Along with the rising of the droplet, both the 
internal and external temperature fields around the droplet surface are re- 
distributed due to the action of heat convection. In the process, the isotherm 
T — moves up and approaches the top of the droplet. At i = 20, a small 
cap-type isotherm T = is formed within the droplet. Meanwhile, under 
the action of heat conduction, the thermal energy is transferred from out- 
side of the droplet to inside, so the temperature inside the droplet increases. 
At t — 60, the temperature of cap-type isotherm in the droplet reaches 
T — 2. Hence, both the temperature fields inside and outside the droplet 
increase with time. Only when both the internal and external temperature 
fields satisfy linear relations with the same slope, does the droplet migration 
reach a steady state. For large Ma numbers, the effect of heat convection is 
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stronger than that of heat conduction. Fig. 9 displays the time evolution 
of isotherms at large Ma (Ma=402.5, Rc=5.93), which corresponds to that 
of migration velocity shown in Fig. 4. In the initial (accelerating) stage of 
migration, the evolution of the isotherms is similar to that shown in Fig. 8. 
In the following (oscillation) stage of migration, the bent isotherms in the 
droplet move to approach the top of the droplet and are converted to pea- 
type ones, and then to earphone-type ones. In the last (accelerating) stage 
of migration, the earphone-type isotherms are transformed into those with 
two symmetric vortices. In the whole migration process, although the tem- 
perature inside the droplet increases, the temperature of minimal isotherm 
is still kept at T = 0. It imphes that the thermal energy transfer from out- 
side of the droplet to inside is weaker than that shown in Fig. 8. Fig. 10 
displays the time evolution of isotherms at a still larger Ma (Ma=2191.6, 
Re=32.3), which corresponds to that of migration velocity shown in Fig. 4. 
In the whole (accelerating) process of migration, the initial evolution of the 
isotherms is similar to that shown in Fig. 9, except for the slower movement 
of the isotherms in approaching the top of the droplet. Then, the pea-type 
isotherms are converted to earphone-type ones. And finally, the isotherms 
with two vortices symmetric about the vertical diameter in the droplet are 
formed. At t — 60, the fact that the minimal isotherm T — has larger 
closed area in the droplet means that the thermal energy transfer from out- 
side of the droplet to inside is weaker than that shown in Fig. 9. Thus, at 
large Ma numbers, although the temperature outside the droplet increases 
fast as the droplet rises, the temperature inside the droplet has only a slow 
increase. The droplet migration does not reach a steady state, and is thus a 
unsteady process. 

To further quantitatively depict the steady and unsteady migration pro- 
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cesses, we will investigate the time evolution of temperature fields inside and 
outside the droplet. Fig. 11 displays the temperature at the point {xc,Zc) 
inside the droplet and the point (xc, Zc + 2) outside the droplet in the migra- 
tion processes, where {xc, is the center of the droplet in the laboratory 
coordinate system. It is evident that for the moderate Ma(=44.9) number 
time evolution curves of the temperature at these two points after to = 20 
are approximately linear and parallel. Both the temperature inside and out- 
side the droplet satisfy the linear relation: Tj = Tj(to) + V^oit — to), which 
indicates a steady migration process with the constant velocity = Voq. 
However, for large Ma numbers, although the time evolution curves of the 
temperature inside and outside the droplet after to = 20 are approximately 
linear, but they are not parallel. The slope of the time evolution curve of the 
temperature inside the droplet is smaller than that outside the droplet, so the 
difference of the temperatures at these two points increases as time increases. 
It implies that the terminal droplet migration does not reach a steady state, 
and is thus a unsteady process. Therefore, the advection around the droplet 
is a more significant mechanism for heat transfer across/around the droplet 
at large Ma numbers. 

3.2.2 Flow field with the temperature gradient G=9 K/cm 

In the space experiments with G=9 K/cm[2Hj, Re and Ma lie respectively 
in the ranges of 3.2-89.8 and 148-4103, their specific values depending on 
Ro- To simulate the experimental processes, the physical coefficients in the 
droplet migration processes are determined by changing Ro. The correspon- 
dence of Re, Ma and Ca to Ro is presented in Table III, where Re is in 
a range of moderate values and Ma have both moderate and large num- 
bers. Fig. 12 displays the time evolution of droplet migration velocities 
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for five sets of non-dimensional coefficients. In the present range of Ma, 
the curves of migration velocities versus time are classified into three types 
based on their characters. At Ma=33.5(Re=0.49), the initial migration ve- 
locity increases sharply near t = 2, and then drops to approach a steady 
value. For Ma=301. 9-838. 6(Re=4. 45-12. 4), after the increasing-decreasing 
oscillation process, the droplet migration is in an accelerating state. For 
Ma=1643. 6-271 7.1(Re=24. 2-40.1), the droplet migration velocity increases 
monotonously with time and decreases with increasing Ma. We can thus con- 
clude that in the time frame under investigation the thermocapillary droplet 
migration is steady at moderate Ma numbers, but becomes unsteady at large 
Ma numbers. To further compare with the experimental investigation |28j . we 
take several fixed migration distances and determine the relation between 
instantaneous non-dimensional migration velocity Vz and Ma numbers. The 
numerical and experimental results are plotted in Fig. 13, from which it 
is evident that both the numerical and experimental migration velocities of 
droplet decrease as Ma increases in the range of large Ma numbers. Hence, at 
large Ma numbers, the above simulation results are in qualitative agreement 
with those of experimental investigations. 

4 Theoretical analysis of thermocapillary droplet 
migration at large Ma numbers 

4.1 Quasi-steady state assumption 

In general, the surface tension is a linear decreasing function of the local 
temperature. For a temperature field with its gradient in the z direction, 
the generated surface tension force is a net force along the surface and the 
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droplet starts to move towards the warm side under the action of net force. 
When the net force acting on the droplet at the flow direction is zero, the 
thermocapillary droplet migration reaches a steady process. However, due 
to the variation of physical parameters with the ambient temperature, the 
migration process may not reach any steady state. Only when the migration 
is sufficiently slow that the order of relevant time scale for the transport 
process to generate steady velocity and temperature flelds is smaller than 
that for the droplet to move an appreciable distance, the assumption of 
the quasi-steady state is valid. It means that after experiencing an initial 
unstable migration process, the droplet migration may reach a steady state 
at the time to and the position ro = ^ok, i.e., migrating with a constant 
droplet migration speed Voo- Using the coordinate transformation from the 
laboratory coordinate system to a coordinate system moving with the droplet 
velocity Ko 

r = f + ro + Ko(^-^o)k, Vi{r,t) ^ Vi{f) + V^k, Ti{r,t) ^ fi{r) + zo + V^{t 

(25) 

the energy equation in Eqs.(6) can be formulated as 

Ko + V • (v,f,) = i^fV • (km), (26) 

in a polar coordinate system f = (^,0). Under the assumption of non- 
deformable droplet, the radial coordinate axis is the outer normal vector of 
the interface. By using the transformation (25), the boundary conditions (7) 
and (8) can be respectively written as follows 

vi = -l/ook,fi^fcose, (27) 



20 



at places far away from the droplet and 



Ti(i,e)-- 



--U2n{l,0) = 0, 

V2s{i,e), 

f2{l,9), 



(28) 



at the interface of the two fluids. Thus, once the droplet migration reaches 
a steady state, the above problem of (6) (7) (8) in the laboratory coordinate 
system can be described by steady energy equations (26) with boundary 
conditions (27) (28) in the coordinate system moving with the droplet velocity. 
This implies the overall steady-state energy balance with two phases in the 
flow domain in the co-moving frame of reference. 

4.2 Nonconservative integral thermal flux across the 
droplet surface at large Ma numbers 

To analyze the energy equations with a small parameter e = l/-\/T4o^Ma, 
Eqs. (26) are rewritten as 



where both the velocity fields Vj in the surrounding fluids and in the droplet 
are rescaled by Voo- 

To confirm the overall steady-state energy balance of two phases in the 
flow domain with respect to the co-moving frame of reference, we have to 
integrate Eqs. (29) (30) with an asymptotic expansion of the outer tempera- 
ture field at infinity with respect to the small parameter e. To determine the 
asymptotic behavior of fi at f 1, we rewrite Eq. (29) as follows 



l + V-(vifi) = e2V-(Vfi), 
l + V-(V2f2) = e'«2V-(Vf2), 



(29) 



(30) 



dTi vi m 



(31) 



dr r do 
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Let 

ui = u1 + 0(1), 

t)i = t)? + o(l), (32) 
fi = fo + o(l), 
we have the energy equation to leading order 

By using the characteristic hne method, the primary approximation of the 
outer temperature field in the continuous phase is derived as 

Ti = fcos^+ \\v\^me -u\Q,o^e -\)lu%df ^o{\), (34) 

where \E' ~ (f — 1/f) sin 6*. For moderate Re numbers, the velocity fields in 
Eq. (34) can be described by the potential fiows. By using the scaled inviscid 
velocity field in the continuous phase fiow passing a circular cylinder 

u° = -cos^(l - ^), 
t;? = sin^(l + ^). 



Eq. (34) is written as 



1 2^V(f - l/f)2 - 1 
r2 - 1 ±^1 - ^2/(f- l/f)2 ' 



Ti = fcos9+ I , _ / ' J ^df + o{l), (36) 



where \l/[= sin0(f— 1/f)] is streamfunction of the continuous phase, the sym- 
bol " +" before the integral is so determined as to preserve the monotonously 
increasing trend of Ti(f, 0) with f(> 1) in the continuous phase and the sym- 
bol " ±" in the integral depends on the value of 6 (the symbol " -|- " /" — " 
corresponds to 6* G [0, 7r/2)/[7r/2, tt)). At f ^ 1, Eq. (36) can be expressed 
as 

= fcos9 — 4 cos 6' + o{l), 
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where ^ ~ sin 9f. 

In physics, the thcrmocapillary migration of a planar droplet has the 
mirror symmetry about the coordinate axis = or tt, so the overall energy 
balance in the whole flow domain e [0, 27r) can be generated by combining 
the energy balance in two connected flow domains 9 e [0, tt) and 9 e [0, tt) 
through the transformation 9 — 2t: — 9. Integrating Eq. (29) and Eq. (30) 
in the continuous phase domain (f G [l,foo],^^ G [0,27r]) and within the 
droplet region (f G [0, 1],^ G [0,27r]), and then transforming them to linear 
integration on the droplet surface and the surface at infinity by using the 
Green's formula, we have 



l) + y" ■UinTi\f^ds- j UinTi\ids = e^{j ^1*^=°^*"/ ^'^^^-^ ^^^^ 
and 



Tx + j> U2nT2\ids = e^K2 j> ^^\ids. (39) 

Using the normal velocity boundary condition at the interface in (28) and 
the temperature field at the infinity in (37), we can derive 

9fi , , TT 1 , , 1 , TT 



dn ^ 



^- = -J(l-^)+^(^)--5 (40) 



and 



-^\^ds = —^. 41 
on K2e 



J an K2e 

To analyze the thermal flux near the boundary, we write the integrals of Eq. 
(40) and Eq. (41) in their discretization forms and simplify the expressions 
in terms of the mirror symmetrical relationship ^|i,6> = ^|i,27r-6> as 

1=1 1=1 



and 



1=1 1=1 
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where — 27r/N. And thus we arrive at a conclusion that there must be 

some interface points 9i G [0,7r] where the following equation holds 

r)T riT 

^(1,^.)<0<^(1,^.) (44) 
or some interface points 6i and Oj G [0, tt] where the following equations hold 

0< #(1,^0 < #(1,^0, 



(45) 

^(1,^,)<^(1,^,)<0. 
Physically, this means that near these points 9i the thermal energy is trans- 
ferred from the interface to outside (the surrounding fluid) as well as from 
the interface to inside (the droplet) or near these points 0i/9j the transfer- 
ence of thermal energy from outside/the interface to the interface/outside is 
weaker /stronger than that from the interface/inside to inside/the interface. 
On the one hand, if Eq. (44) can satisfy the thermal flux boundary condi- 
tion in Eqs. (28), thermal sources inside the interface will be introduced to 
balance the transference of thermal energy. On the other hand, if Eqs. (45) 
can satisfy the thermal flux boundary condition in Eqs. (28), thermal sinks 
inside the interface or thermal sources in the droplet will be introduced to 
decrease the transference of thermal energy from the interface to outside or 
increase the transference of thermal energy from inside to the interface. Since 
there is absolutely no thermal sources or sinks inside the interface or ther- 
mal sources in the droplet, the above transport processes of thermal energy 
near the interface seem impossible. It means that the thermal flux across the 
droplet surface is nonconservative. Moreover, from Eq. (40) and Eq. (41), 
we have 

[A:2^|i - ^^Ws = ^(1 + ^) = n{l + '^)V^Ma. (46) 
or or K2 K2 

Since both k2 and K2 are positive, we have 

k2 f ^\ids > f ^\ids. (47) 



dr J dr 
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From Eqs. (28), we obtain the equivalent integral thermal flux at the bound- 
ary 



So, if the thermal flux boundary condition in Eqs. (28) is satisfled, Eq. (47) 
should be reduced to Eq. (48), which seems impossible. It is termed as 
a nonconservative integral thermal flux across the surface for the steady 
thermocapillary droplet migration at large Ma numbers. This implies the 
overall steady-state energy unbalance of two phases in the flow domain in 
the co-moving frame of reference. 

Eq. (47) indicates that at large Ma numbers the integral thermal flux 
across the surface within the droplet is larger than the surface thermal flux 
with respect to the continuous phase fluid. However, it should be noted that 
the droplet migration is, at that time, still in an unsteady state. In the 
analytical and numerical results [26| |27] . the steady migration velocity (in 
direct proportion to Ma) is large at high Ma numbers. Under the condition 
of large migration velocity, it is unlikely that the order of relevant time scale 
for the transport process to generate steady velocity and temperature flelds 
is smaller than that for the droplet to move an appreciable distance. So, 
due to the variation of physical parameters with the ambient temperature, 
a steady migration process may not be reached. Both experimental results 
in Fig. 4(a) (b) of [II] and Fig. 4(b) (d) of [28] clearly display that the 
thermocapillary droplet migration at large Ma numbers is in an accelerating 
state and does not reach any steady one. Moreover, numerical simulations 
of the thermocapillary motion of deformable and non-deformable droplets in 
[I8j and the above section 3 indicate that the assumption of quasi-steady 
state is not valid for large Ma numbers. Thus, it is clear that the invalid 
assumption of quasi-steady state for the thermocapillary droplet migration 




(48) 
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process is a reasonable explanation for the nonconservative integral thermal 
flux across the droplet surface. 

5 Conclusion and discussions 

In this paper, numerical studies are carried out for thermocapillary migra- 
tion of a planar non-deformable droplet in two uniform temperature gradi- 
ents at moderate and large Ma numbers by using the front-tracking method. 
Some calculations at moderate and large Ma numbers are performed to an- 
alyze the thermocapillary migration for droplets with different sizes. In the 
range of droplet radius under study, Re takes moderate values, and thus the 
computed flow flelds are similar. There appear different types of migration 
processes in the time frame under investigation depending on the values of 
Ma, which varies within a large range. At moderate Ma numbers, after an 
increase-decrease process in the time evolution of droplet velocity, the droplet 
migration reaches a steady state. In the range of large Ma numbers, the oscil- 
lation process in the time evolution of droplet velocity is transformed into a 
monotonous accelerating process as Ma increases. The terminal droplet mi- 
gration is in an acceleration process and doesn't reach any steady state. The 
instantaneous migration velocity at a flxed migration distance decreases with 
increasing Ma number. The numerical simulation results are in qualitative 
agreement with experimental ones. 

Moreover, in comparing the variations of temperature fields inside and 
outside the droplet, it is evident that at large Ma numbers the weak transport 
of thermal energy from outside of the droplet into inside cannot meet the 
requirement put forward by the steady migration process, which implies that 
the advection around the droplet is a more significant mechanism for heat 
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transfer across/around the droplet at large Ma numbers. 

Furthermore, from the condition of overall steady-state energy balance in 
the flow domain, we have identified a nonconservative integral thermal flux 
across the surface for a steady thermocapillary drop migration in a uniform 
temperature gradient at large Ma (Re) numbers. It may well result from the 
invalid assumption of quasi-steady state, and this conclusion implies that the 
thermocapillary drop migration at large Ma (Re) cannot reach any steady 
state and is thus a unsteady process. 

We emphasize that all of the above numerical and theoretical results 
about the thermocapillary migration system of droplets involves assump- 
tions of planar non-deformable interfaces and is subject to constant physical 
parameters. As mentioned in the section 1, the simple modelling is easily 
applied to explore the dynamical mechanism of droplets, but has potential 
drawbacks to reproduce the experimental results |28]. Extension to the more 
realistic case of three-dimensional deformable droplets migrating in a flow 
field with temperature-dependance physical parameters remains to be imple- 
mented. 
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Table I. Physical parameters of the continuous fluid (5cst Silicone oil) 
and the droplet (FC-75) at temperature 25°C, which are the working media 



in the space experiment pH]. 




p{g/cm^) 


/i(10 "^dyns/crn^) 


k{W/mK) 


k(10 '^crn^/s) 


Silicone oil 


0.91 


4.268 


0.111 


6.915 


FC-75 


1.77 


1.416 


0.063 


2.018 



Table II. Correspondence of non-dimensional parameters Re, Ma and Ca 
to the droplet radius Rq for the droplet migration in a flow fleld with the 
temperature gradient G=12K/cm. 



-Ro(cm) 


Re 


Ma 


Ca 


0.05 


0.66 


44.7 


0.0044 


0.15 


5.93 


402.5 


0.013 


0.25 


16.5 


1118.1 


0.022 


0.35 


32.3 


2191.6 


0.031 


0.45 


53.4 


3622.8 


0.040 



Table III. Correspondence of non-dimensional parameters Re, Ma and Ca 
to the droplet radius Rq for the droplet migration in a flow fleld with the 
temperature gradient G'=9K/cm. 



-Ro(cm) 


Re 


Ma 


Ca 


0.05 


0.49 


33.5 


0.0033 


0.15 


4.45 


301.9 


0.010 


0.25 


12.4 


838.6 


0.017 


0.35 


24.2 


1643.6 


0.023 


0.45 


40.1 


2717.1 


0.030 
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Figure caption 

Fig. 1. Schematic of tlie computation domain for a planar droplet mi- 
gration. The top and bottom walls are non-slip boundaries and the left and 
right boundaries are periodic ones. 

Fig. 2. Droplet migration velocity versus non-dimensional time for three 
grid resolutions 64 x 128, 96 x 192 and 128 x 256 at a fixed domain 4x8 
under Re=5, Ma=20, Ca=0. 01666 and pij p\= ^2/^1= ^2/^1= /t2//«i=0.5. 

Fig. 3. Time evolution of the droplet migration velocity for grid resolution 
96 X 192 and its comparison with Nas & Tryggvason's result [15] for the same 
parameters as given in Fig. 2. 

Fig. 4. Droplet migration velocity in a fiow field with the temperature gra- 
dient G = 12K/cm versus non-dimensional time at Ma=44.7, 402.5, 1118.1, 
2191.6 and 3622.8. 

Fig. 5. Instantaneous thermocapillary migration velocity of the droplet in 
a fiow field with the temperature gradient G = 12K/cm at a fixed migration 
distance = l/1.5/2cm denoted by diamonds/deltas/squares versus large 
Ma numbers. The experimental results {28] rescaled by Vygb/vq = 2/[(2 + 
3/i2//ii)(2 + k2/ki)] are plotted and denoted by circles. 

Fig. 6. Computed velocity fields at t=20 under i?o=0.25cm, Re=16.5, 
Ma=1118.1 in (a) the laboratory coordinate frame and (b) the reference 
frame moving with the droplet. 

Fig. 7. Streamlines in a reference frame moving with the droplet under 
i?o=0.25cm, Re=16.5, Ma=1118.1. Their time evolution is displayed in 5 
small figures from left to right. The non-dimensional time is chosen as 3, 10, 
20, 40 and 60, respectively. 

Fig. 8. Isotherms in a laboratory coordinate frame are selected from the 
computation of the droplet migration under i?o=0.05cm, Re=0.66, Ma=44.7. 
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Notation is the same as in Fig. 7. 

Fig. 9. Same as Fig. 8, except i?o=0.15cm, Re=5.93, Ma=402.5. 

Fig. 10. Same as Fig. 8, except /2o=0.35cm, Re=32.3, Ma=2191.6. 

Fig. 11. Time evolution of temperature at point (xc, Zc) inside the droplet 
and point (xc, Zc + 2) outside the droplet, where (xc, Zc) is the center of the 
droplet in the laboratory coordinate system. 

Fig. 12. Droplet migration velocity in a flow field with the temperature 
gradient G = 9K/cm versus non-dimensional time at Ma=44.7, 301.9, 838.6, 
1643.7 and 2717.1. 

Fig. 13. Instantaneous thermocapillary migration velocity of the droplet 
in a flow field with the temperature gradient G = 9K/cm at a fixed migration 
distance = 1/1. 5 /2cm denoted by diamonds/deltas/squares versus large 
Ma numbers. The experimental results |28j rescaled by Vygb/vq = 2/[(2 + 
3/i2//^i)(2 + k2/ki)] are plotted and denoted by circles. 
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